tic

%%% This file is to put different theta cases together and plot graphs for easy comparison.

clear all
close all

global parms i kr Hr jr cr Sf kf Nf cf  sigmaf T0 k0


n0 = 0.00832913;      % Initial level of mining investment
c0 = 0.6619974; % Initial value of consumption for calculate lambda_0
k0 = 3.6071282734; % Initial value of K for the differential equation

% ######################
% ## No externality   ##
% ######################
% Tbar = 366;
% k_Tbar = 133937.035084899981;
% S_T0 = 1595.42705305;
% T0 = 98, k0 = 3.607023,N0 = 9.2457e-05, S0 = 20.4151
% 
% 
% #######################################
% ##  rho = 0.05*psi, theta = 0.95*psi ##
% #######################################
% Tbar = 387;
% k_Tbar = 136086.903981801035;
% S_T0 = 1582.48792067;
% T0 = 103, k0 = 3.608044,N0 = 8.9564e-05, S0 = 20.4153
% 
% 
% #######################################
% ##  rho = 0.25*psi, theta = 0.75*psi ##
% #######################################
% Tbar = 400;
% k_Tbar = 146197.463759450009;
% S_T0 = 1579.61845203;
% T0 = 106, k0 = 3.607139,N0 = 3.2377e-05, S0 = 20.4068
% 
% #######################################
% ##  rho = 0.95*psi, theta = 0.05*psi ##
% #######################################
% Tbar = 436;
% k_Tbar = 204550.525554385968;
% S_T0 = 1570.19594116;
% T0 = 112, k0 = 3.601186,N0 = 9.5647e-05, S0 = 20.4126


% theta = [no externality; thetaj=0.05psi; thetaj = 0.25psi; thetaj = 0.95psi]
theta = [0, 0.05*0.33, 0.25*0.33, 0.95*0.33; 0, 0.95*0.33, 0.75*0.33, 0.05*0.33];
Tbar = [366, 387, 400, 436]; 
k_Tbar = [ 133937.035084899981, 136086.903981801035, 146197.463759450009, 204550.525554385968];
S_T0 = [1595.42705305,1582.48792067, 1579.61845203, 1570.19594116];


% theta = [no externality; thetaj = thetaB = 0.5; thetaj = 0.25; thetaB = 0.05]
% theta = [0, 0.5*0.33, 0.25*0.33, 0.95*0.33; 0, 0.5*0.33, 0.75*0.33, 0.05*0.33];
% Tbar = [394, 421, 400, 436]; 
% k_Tbar = [ 133937.034433349967,  161864.0039534, 146197.463759450009, 204550.525554385968];
% S_T0 = [1579.33947006,1569.4260983, 1579.61845203, 1570.19594116];
scenario = [theta;Tbar;k_Tbar;S_T0];
[Tmat0, wholemat0, renewmat0, fossilmat0] = mainfun(scenario(:,1));
[Tmat1, wholemat1, renewmat1, fossilmat1] = mainfun(scenario(:,2)); 
[Tmat2, wholemat2, renewmat2, fossilmat2] = mainfun(scenario(:,3));
[Tmat3, wholemat3, renewmat3, fossilmat3] = mainfun(scenario(:,4)); 
save results.mat